Under consideration for publication in J. Fluid Mech. 



1 



Transition in pipe flow: the saddle structure 
on the boundary of turbulence 

Y. DUGUET,H A.P. WILLIS M 
AND R.R. KERSWELLi 

^ School of Mathematics, University of Bristol, BS8 ITW Bristol, United Kingdom 

(Received November 13th 2007 and in revised form May 27th 2008) 

The laminar-turbulent boundary E is the set separating initial conditions which re- 
laminarise uneventfully from those which become turbulent. Phase space trajectories on 
this hypersurface in cylindrical pipe flow look to be chaotic and show recurring evi- 
dence of coherent structures. A general numerical technique is developed for recognising 
approaches to these structures and then for identifying the exact coherent solutions 
themselves. Numerical evidence is presented which suggests that trajectories on E are 
organised around only a few travelling waves and their heteroclinic connections. If the 
flow is suitably constrained to a subspace with a discrete rotational symmetry, it is pos- 
sible to find locally-attracting travelling waves embedded within E. Four new types of 
travelling waves were found using this approach. 



1. Introduction 

Transition to turbulence in cylindrical pipe flow is governed by one single dimensionless 
parameter, the Reynolds number Re := UD/v where U is the mean flow speed along the 
pipe, D the pipe diameter and v the kinematic viscosity of the fluid (Reynolds 1883). 
Despite the simplicity of the set-up, the reason for transition remains obscure due to the 
linear stability of the laminar Hagen-Poiseuille flow (Hagen 1839, Poiseuille 1840) and 
the sensitivity of the process to the exact shape and amplitude of disturbances present. 
In most experiments, transition is observed at Re ~ 2000 (e.g. Wygnanski & Champagne 
1973) but can be triggered as low as Re — 1750 (Peixinho & MuUin 2006) or delayed to 
Re = 100,000 in very carefully controlled experiments (Pfenniger 1961). Until recently, 
the only firm theoretical result was the energy stability bound of Re = 81.49 (Joseph & 
Carmi, 1969) below which all disturbances are guaranteed to decay monotonically. This 
is, however, more than an order of magnitude below the observed value for transition. 

An important step forward in understanding the transition process was the discovery of 
disconnected solutions to the Navier-Stokes equations in a cylindrical pipe (Faisst & Eck- 
hardt 2003, Wedin & Kerswell 2004, Kerswell 2005, Pringle & Kerswell 2007). These exact 
solutions are travelling waves (TWs) which appear through saddle-node bifurcations as 
in other shear flows (Nagata 1990, Waleffe 1997,1998,2001,2003). All these solutions are 
linearly unstable though they have a very low-dimcnsional unstable manifold. There is 
much interest in these solutions, as very similar structures have been observed transiently 
in experiments (Hof et al. 2004, 2005) and direct numerical simulations (Kerswell & Tutty 
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2007, Schneider et al. 2007a, Willis & Kerswell 2008a). These states can generally be di- 
vided into 'upper-branch' and 'lower-branch' TWs, based on whether they have high or 
low wall shear stress. Lower branch solutions are believed to sit on a hypersurface that 
divides phase space into two regions: one where initial points lead directly to the lami- 
nar state, the other where initial conditions lead to turbulent episodes (Kawahara 2005, 
Wang et. al. 2007, Kerswell & Tutty 2007, Viswanath 2008). The simple translational 
behaviour of these travelling wave solutions is inherent to the method used to find them, 
and undoubtedly masks an even larger variety of more complex exact solutions. 

The boundary between laminar and turbulent trajectories - labelled S hereafter - is 
formally a separatrix if the turbulent state is an attractor. At low Re, however, turbu- 
lence may ultimately decay after a long transient in which case the laminar state is the 
unique global attractor. Then the boundary S is generalised to the; dividing set in phase 
space between trajectories which smoothly relaminarise and those which do undergo a 
turbulent evolution. E is thought to be of codimension 1 in phase space but one can 
a priori not exc;Iude a more complex fractal structure, as suggested by the dependence 
of lifetime on initial conditions in simulations of pipe flow (Faisst & Eckhardt 2004). 
Trajectories which start in S - hereafter called the 'edge' following Skufca et al. (2006) - 
stay in E for later times by definition and hence the long time dynamics are of obvious 
interest. The long time behaviour on E has already been found to be a periodic solution 
in plane Poiseuille flow in a pioneering study by Toh & Itano (1999, see also Itano & 
Toh 2001). Skufca et al. (2006) studied a 9-dimensional model of plane Couette flow 
(PCF) to reveal an attracting periodic orbit at low Re and an apparently chaotic state 
at higher Re. However, recent fully-resolved simulations have shown that the asymptotic 
behaviour is an attracting TW in PCF (Schneider et al. 2008, Viswanath 2008) and a 
chaotic attractor in a short cylindrical pipe of length L = 5D (Schneider et al. 2007b). 
Interestingly, this seemingly chaotic end state looks to be centered around an 'asymmet- 
ric' TW (Pringle & KersweU 2007, MeUibovsky & Meseguer 2007). 

The purpose of this paper is to explore the dividing hypersurface E in pipe flow with 
the following objectives: 

(a) to establish that the dynamics restricted to this laminar-turbulent boundary ex- 
plores many different saddle points embedded in it; 

(6) to find evidence for heteroclinic or 'relative' homoclinic connections between these 
saddle points; 

(c) to explore E restricted by a discrete rotational symmetry in order to ascertain 
whether the limiting behaviour remains chaotic or can be a simple attractor; 

(d) to develop a practical and general way to find TWs and periodic orbits without 
detailed knowledge of their spatial structure. 

The works cited above employ a reduced computational domain for their calculations, 
imposing short-wavelength periodicity in one direction (along the pipe, Schneider et al. 
2007b) or two (in the spanwise and streamwise directions in plane Couette flow, Schneider 
et al. 2008, Viswanath 2008). This is very much in the spirit of the Minimal Flow Unit 
pioneered by Jimenez & Moin (1991) and used with considerable success to identify key 
mechanisms underpinning turbulence (Hamilton et al. 1995). In this paper, we also adopt 
this same approach by concentrating on pipes up to 5 D long. The ensuing reduction in 
the degrees of freedom of the flow allows a much more detailed exploration of the flow 
dynamics albeit restricted to a siibspace of strict spatial periodicity. This is a necessary 
preliminary to motivate a more carefully focussed study in much longer pipes which can 
capture spatially-localised turbulent structures (e.g. Priymak & Miyazaki 2004, Willis & 
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Kerswell 2007, 2008a) but at considerable numerical cost still. Flow structures found in 
the short pipes considered here also exist, of course, in a longer pipe but will be more 
unstable due to the greater degrees of freedom present there. It is also worth remarking 
that a recent study (Nikitin 2007) has shown how pipe flow is slow to lose its spatial 
periodicity when this is restriction is lifted. 

This paper is organised as follows. Section 2 discusses the formulation and numerical 
methods used to simulate the flow in a pipe and to extract exact recurrent flow solutions. 
Section 3 presents the results obtained using these in 7 subsections. Subsection 3.1 recalls 
the method used to follow trajectories on S and confirms that the limiting behaviour 
looks chaotic (Schneider et al. 2007b). Sections 3.2 and 3.3 discuss how near-recurrent 
states are identified and using a Newton-Krylov algorithm demonstrate that a number of 
unstable travelling waves are approached. Subsection 3.4 offers evidence for the existence 
of a 'relative' homoclinic connection between a travelling wave and the same wave rotated. 
Subsection 3.5 looks for recurrent flow structures in S restricted by a discrete rotational 
symmetry and identifies new exact travelling wave solutions. Subsection 3.6 shows that 
within this subspace the limiting state of S is a simple attractor. Subsection 3.7 confirms 
that all the TWs found in the earlier subsections are indeed embedded in the laminar- 
turbulent boundary. The paper ends with a discussion in section 4. 



2. Numerical Procedure 

2.1. Governing Equations 

We consider the incompressible flow of Newtonian fluid in a cylindrical pipe and adopt 
the usual set of cylindrical coordinates (s, 9, z) and velocity components u = us+vd+wz. 
The domain considered here is {s,9,z) E [0:1] x [0:27r] x [0:i], where L = 27r/a is the 
length of the pipe and lengths are in units of radii {D/2). The flow is described by the 
incompressible three-dimensional Navier-Stokes equations 

-^ + (w V)M = -Vp+— V^u, (2.1) 

V-M = 0. (2.2) 

The flow is driven by a constant mass-fltix condition, as in recent experiments (e.g. 
Peixinho & Mullin 2006). The boundary conditions are periodicity across the pipe length 
u{s, 6, z) = u{s, 6,z + L) and no-slip on the walls u{l, 6, z) = 0. 

2.2. Time-stepping code 

The basic tool for the numerical determination of exact recurrent states is the accurate 
time-stepping code used by Willis & Kerswell (2007,2008a). The velocity fleld is derived 
from two scalar potentials ^ and 

w = V X (^-l) + V X V X ($1) (2.3) 

and the incompressible Navier-Stokes equations are re-written using the formulation in- 
troduced by Marques (1990). The two scalar potentials are discretised using high-order 
finite differences in the radial direction s and spectral Fourier expansions in the azimuthal 
direction and axial direction z. For example, the decomposition of the scalar potential 
$ at a radial location Sj, {j = 1, ...,N), reads 

K M 

^Sj,e,z,t;a,mo)= ^ ^ $,.,„(t)e'(™o'»«+"'=-) (2.4) 

k=—K m= — M 



4 



Y. Duguet, A. P. Willis and R.R. Ker swell 



The positive integer toq refers to the discrete rotational symmetry 

27r 

R„o : {u,v,w,p){s,9,z) {u,v,w,p){s,9 -\ , z) (2.5) 

mo 

of the flow (trivially toq = 1 means no rotational symmetry is imposed). The resolution 
of a given calculation is described by a vector {N, M, K). It is adjusted until the energy 
in the spectral coefficients drops by at least 5 and usually 6 decades from lowest to 
highest-order modes. The time stepping is 2nd order accurate with At updated using an 
adaptive method based on a CFL condition. 

The number of (real) degrees of freedom, which defines the dimension of phase-space, 
is 0{8MNK). The set of all complex coefficients X = {^jkm, '^jkm} defines our phase- 
space with its usual Euclidean norm |X| = V X ■ X . Symmetries exhibited by TWs such 
as the shift-and-reflect symmetry 

7T 

S : {u,v, w,p){s, 9, z) {u, —v, w,p){s,—9, z -{ ), (2-6) 

a 

and the mirror symmetry (e.g. Pringle & Kerswell 2007, Pringle et al. 2008) 

Z : {u,v,w,p){s,9,z) {u,-v,w,p){s,-9,z), (2.7) 
are not imposed in the simulations. 

2.3. The Newton- Krylov method 

The spectral expansion defined above converts the Navier-Stokes equations into an au- 
tonomous dynamical system of the form 

^ = F(X). (2.8) 

Travelling wave solutions are 'relative equilibrium' solutions of the Navier-Stokes, steady 
in an appropriate Galilean frame. It is numerically convenient to seek and interpret them 
as a special case of "relative periodic orbit" (RPO) (Viswanath 2007, 2008). We hence 
developed an algorithm to look for RPOs defined as zeros of the functional 

g^\X{T)-'^'--'^'^ -X{0)\\ (2.9) 

Here, X(T) is the point at time T on the trajectory starting at time t = from the point 
X{0). x~^^'~^^ is the point in phase space corresponding to the state X shifted back 
in space by the distance Az in the axial direction and by the angle A^ in the azimuthal 
direction. A shift by (Az, AS) corresponds in phase-space to the transformation 

(*jfe™, <i>jfe™)e-'(""°^^+"'=^^) ^ (^,kra,<^,kra). (2.10) 

A zero of g corresponds to a flow repeating itself exactly after a time T, but at a different 
location defined by Az and AS. A travelling wave solution is a RPO where there is a 
degeneracy between the shifts Az and AS, and the period T. For the case of a TW 
propagating axially with speed c, Az — cT, where T is an apparent period. To remove 
this degeneracy, we impose Az = L unless otherwise specified. As the majority of known 
TWs do not rotate, AS = is also assumed (those that do rotate do so rather slowly, 
Pringle & Kerswell, 2007). Starting from a good initial guess for X{0) and an estimate 
of the period T (to be discussed in Section l3.2p . we minimise the residual g using a 
Newton- Krylov algorithm, based on a GMRES algorithm (Saad & Schultz, 1986). The 
size of the dynamical system (typically 10^ degrees of freedom) necessitates the use of a 
matrix-free formulation. The use of an inexact Krylov solver also leads to a significant 
gain in computation time (we used the choice 2 in Eisenstat & Walker, 1995). Moreover, 
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we embed the Newton solver into a more globally convergent strategy in order to improve 
likeliness of convergence, using a double dogleg step technique (Dennis & Schnabel, 1995, 
Brown & Saad 1990). For the special case of TWs with a short period, we assume full 
convergence when the normalised residual ^/\X{Q)\ is less than O(10~^"). 

2.4. Stability of Travelling Waves 

Once a travelling wave solution is known, along with its axial propagation speed c, it 
can be expressed as a steady solution in the frame moving at speed c. In this Galilean 
frame, the stability of the solutions can be studied numerically using an eigenvalue solver 
based on an Arnoldi algorithm. This yields the leading eigenvalues whose real part, when 
positive, indicates the growth rate of infinitesimal perturbations to the exact solution in 
the moving frame. 

3. Results 

We present 4 'edge' calculations, each motivated by a different question. In the first 
(see ^3.H -3). the edge in a pipe length L « 5D (a = 0.625) is calculated starting from 
a typical turbulent initial condition in order to investigate whether there are coherent 
structures buried in it. We are able to confirm the presence of what looks to be a chaotic 
attractor as found recently by Schneider et. al. (2007b). In the second edge calculation 
( §3.4p we look for numerical evidence of heteroclinic connections between saddle points by 
using a perturbed TW as an initial condition. In the third ( ^3.5p we impose the discrete 
rotational symmetry R2 on the flow to exploit the saddle structure of the subset of S to 
reveal new TWs which possess R2-symmetry. In the fourth fi j3.6[) we look for evidence 
of multiple attracting TWs to demonstrate that under a rotational symmetry constraint 
the large-time dynamics need not be chaotic. 

3.1. Calculating Edge Trajectories 

In this subsection we choose a pipe with a — 0.625 {L = 5.02651?) and Re = 2875, set 
mo = 1 so that there is no restriction on the rotational symmetry, and take a numerical 
resolution of (30,15,15). For this value of Re, turbulence, when triggered, is sustained 
over the typical simulation time which is much longer that the time taken to relaminarise. 
In order to constrain a numerical trajectory to stay on the edge surface, S, we use a 
shooting method analogous to that first used in plane Couette flow by Toh & Itano 
(1999). We first produce a long turbulent trajectory and pick a state u* of relatively 
low three-dimensional disturbance energy E^dit) (defined as the kinetic energy in the 
streamwise-dependent Fourier modes fc 7^ 0). Initial states are then chosen from the line, 

up:^{u*)+f3{u*-{u*)), (3.1) 

parameterised by the real number (3 and ( ) :— ^ J^" {..)dz is an axial average. The 
initial condition tt/3=o is two-dimensional and hence does not trigger turbulence, whereas 
tX/3=i is the original turbulent state. This /3-interval is then repeatedly halved with, at 
each stage, the new interval selected such that an initial condition corresponding to /3 
at the lower limit relaminarises, whereas that at the upper limit leads to large turbulent 
values of E^d- /3 is refined close to machine accuracy using double precision arithmetic, 
forcing the energy of the trajectory to stay at an intermediate level corresponding to 
the laminar-turbulent boundary for typical times of O{2Q0 D/U). At the values of Re 
used here, the energy of this boundary is clearly distinct from those associated with 
turbulence and thus makes determination of p unambiguous. Once machine precision 
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has been reached, the process is restarted from nearby states towards the end of the 
trajectory originating from neighbouring [3 values. 

The fuU trajectory, of duration 0(500 D/JJ), does not show any sign of convergence 
towards a simple state but rather displays erratic and unpredictable dynamics, as found 
by Schneider et al (2007b). We deliberately avoid the word 'turbulent' as this refers to 
larger values of Ej,d. Figure [T] shows the time evolution of the E^d{t)- 



3.2. Near- Recurrences on an Edge Trajectory 

Despite the lack of regularity of the energy signals, inspection of all velocity components 
at several random locations in the pipe indicate clearly that the flow on E is sometimes 
nearly periodic in time on short intervals. For the case (a, i?e,TOo) — (0.625,2875, 1), it 
is possible by eye to identify several temporal windows where all velocity components 
approximately oscillate with a given frequency (see Figure [2]) . Within these windows the 
dynamics on the edge appear to be temporally correlated with a short period. A large 
number of snapshots were taken at times = 1,2,3... across the whole trajectory. 
To investigate the possibility of recurrence in the flow, each snapshot state X{ti) was 
used as an initial condition for time-stepping. The normalised distance in phase space 
between this initial point and its evolution in time was then examined for local minima. 
Specifically, we define the scalar residual function 

\X-^'-^0(t)~ X{t,)\ , , 

where Az is a distance by which the state X{t) is shifted back in z for comparison. We 
chose Az = L so that a value of ri (t) = for some t > ti means that the flow in the pipe 
is exactly the flow at time ti. In this case, X{ti) lies on a periodic orbit of the system 
with the time interval t — ti being a multiple of the period (for Az 7^ 0, a vanishing 
residual would indicate a relative periodic orbit). For the ease of calculation, most of 
the time we chose to neglect the possibility of recurrence occurring after a shift in the 
azimuthal direction 9, and therefore A9 = was imposed. Typical plots of ri{t) starting 
from different values of ti are displayed in Figure [31 Every function attains a smallest 
local minimum which is defined as 

Or ■ 

r„iin{ti) := mm{ri{t) : = 0}. (3.3) 

Figure S] shows how rmin varies with the starting point. Phases for which rmin is 
small (from experience below 0.1) may be interpreted as approximate approaches to pe- 
riodic orbits of the system by the edge trajectory. The alternating pattern of maxima 
and minima in rmin is then the signature of a repellor - a set of states which attract 
trajectories before ultimately repelling them away. For hyperbolic dynamical systems, 
trajectories are attracted towards one of these states along their stable manifold and 
ejected away along their unstable manifold (see the sketch in Figure [5]). For the parame- 
ters here, 6 dips in the r^in-proiile are suggestive of 6 approaches, denoted respectively 
by Al, B1,C1, Dl, El, Fl. Parts of the trajectory linking one state to the next (e.g. 
Al Bl or Dl El) may be located in the vicinity of heteroclinic connections be- 
tween the two states, should such a connection exist, or a relative homoclinic connection 
for two successive states which are the same but for a shift. The notion of 'vicinity' de- 
pends on the choice of the norm in phase space. Here, a pragmatic approach is adopted: 
we use the expression 'X is close to a periodic orbit Y' to mean that the Newton-Krylov 
algorithm converges to the periodic orbit Y starting with the initial guess X. 
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3.3. Exact Coherent Structures in E 

We now analyse the neighbourhood of the points X(ti) yielding a low value of rmin 
looking for exact recurrent states for which rmin exactly vanishes. Starting from such 
initial guesses, the Newton-Krylov algorithm defined in !j2.3l was used to further minimise 
fmin ■ Of the six different starting guesses A to F introduced in the previous subsection, 
excellent convergence was obtained in the cases Al, J51, Dl and Fl, with rmin being 
reduced down to O(10~^^). The converged states are labelled respectively Al_0.625, 
i31_0.625, £'1_0.625 and Fl_0.625 in order to distinguish them from the starting points 
Al, Fl, and to indicate the parameter a. The cases CI and El displayed only partial 
convergence to respectively O{10~^) and O(10~^) and it is not possible to say if there 
really are zeros of rmin in this neighbourhood. Despite 'globalisation' improvements to the 
Newton-Krylov algorithm to improve convergence, sample starting guesses away from the 
candidates A, ...,F failed to converge. The converged states Al_0.625, 51-0.625, -Dl_0.625 
and Fl_0.625 all correspond to travelling wave solutions. Close inspection of their spatial 
structure and dynamics shows that they are, in fact, the same travelling wave solution 
modulo a shift in the azimuthal direction: see Figure [HI This state corresponds exactly to 
the 'asymmetric' TW identified in Pringle & Kerswell (2007) whose ^-averaged velocity 
field is strongly reminiscent of the chaotic state calculated by Schneider et al. (2007b). 
This resemblance has also been noted by Mellibovsky & Meseguer (2007) who used a 
different approach to infer the significance of this same asymmetric TW. Interestingly, 
the present study at these parameter values did not find evidence of states with higher 
rotational symmetry in spite of the fact that such TWs are known to be embedded in 
the edge (see fig. 4 of Kerswell & Tutty, 2007). 

The search for recurrent states described above was initially undertaken with A9 = 
disallowing azimuthal propagation of the recurrent patterns. In the case Al, we experi- 
mented with allowing Ad to be updated by the Newton-Raphson scheme. This involved 
solving an additional equation dg/dAO — (Viswanath 2007). Using Newton-Krylov 
with the same starting guess Al but allowing the shift AO to be updated, produced con- 
vergence to another state A'_0.625 which is distinct from yll_0.625. This state, despite 
a velocity profile very analogous to that of y41_0.625, rotates by an angle AO — —0.32 
degrees after travelling one pipe length. This situation is reminiscent of the bifurcation 
diagram obtained by Pringle & Kerswell (2007) for the asymmetric TW (though here 
Re is much higher): for a given Reynolds number, a branch of solutions with helicity is 
connected to the mirror-symmetric branch, and intersects the non-helical subspace twice. 
The solutions at the crossing points are non-helical but nevertheless possess a rotational 
propagation speed ce ^ 0, like our solution yl'_0.625. The fact that such a solution has 
been found here could be explained by the fact that the numerical code used does not 
allow helicity, hence the Newton-Krylov algorithm has no choice but to look for the in- 
tersection of the helical branch with a non-helical subspace. However, other attempts to 
find exact recurrent patterns with rotation invariably converged back to a non-rotating 
{AO = 0) TW. 

3.4. Search for Heteroclinic Connections 

The results so far indicate the possible existence of heteroclinic trajectories linking the 
exact states found. When the two consecutively visited states are the same modulo a 
shift in the azimuthal direction, "relative" homoclinic connection is a more appropri- 
ate term. In an attempt to find further evidence for such connections, a series of initial 
conditions sampling the unstable manifold of a TW were explored. The asymmetric trav- 
elling wave of Pringle & Kerswell (2007) was chosen (i.e. ^1_0.625) with the parameters 
{a. Re, mo) = (0.625, 2875, 1) and orientated so that it was S-symmetric about the = 
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plane. In the corresponding S-symmetric subspace, the solution has exactly two unstable 
eigendirections, denoted by ei and 62, ei being the most unstable one (there are no un- 
stable directions in the dual space of S-anteymmetric perturbations). These two vectors 
were normalised to have unit energy and an initial condition was defined as 

X(0) = Xtw + e (cos ei +sin(/)e2) . (3.4) 

Here, e is a positive parameter chosen small enough that X{Q) is effectively in the unstable 
manifold of Xtvk (e = 10~^ suffices). The angle </> defines a plane spanned by ei and 62 
and vanishes in the direction of ei. For some values of the angle (j), the trajectory starting 
at X(0) returns to the laminar state whereas for others, it becomes turbulent. Boundary 
values of </> between these two behaviours lead to intermediate dynamics demonstrating 
that Xtw sits on the edge E. In all cases the trajectory starts in the S-symmetric 
subspace, but after a finite time this symmetry is broken by numerical noise. A bisection 
method is used to determine (to machine precision) an angle for which the trajectory 
neither evolves to the laminar or turbulent state (see Figure[7]for a sketch of the method). 
The energy contained in the axially-dependent modes is displayed in Figure [5] and the 
residual function Tmin (t) in Figure [D These show that rmin increases exponentially but 
the energy changes little until t ~ 150D/U. After this linear regime, rmin continues to 
increase, then drops to a local minimum value of 0.09 at i ^ 250D/U before increasing 
again. When the state at the local minimum (B4 in Figure [9]) was used as a starting 
point for the Newton-Krylov algorithm, rmin smoothly decreased down to O(10~^^). 
The exact solution found was the asymmetric travelling wave again but now shifted by 
an angle of 51.56 degrees. On the basis that the Newton-Krylov method usually needs to 
be in the vicinity of an exact solution to converge to it, this suggests the existence of a 
relative homoclinic connection near the computed trajectory. Collecting firmer evidence 
is difficult unless the dynamics is further restricted. Recently, significant progress has 
been made along these lines in plane Couette flow by Gibson it et al. (2008) who have 
demonstrated the existence of a heteroclinic connection between two steady solutions on 
E for a minimal flow unit. 



3.5. E under ^12- Symmetry 

The technique developed above to identify recurrent states in the edge can also equally 
be applied to the dynamics within a symmetry subspace. Restricting the flow to be Rm- 
symmetric, for example, improves the possibility of an edge trajectory coming close to 
some of the Rm-symmetric lower branch TWs found to be in the edge by Kerswell & 
Tutty (2007). There is also, of course, the possibility of discovering unknown branches 
of TWs. Attention was restricted to the R2-symmetric subspace by setting mo = 2 in 
the flow representation (|2.4p . At Re = 2875, a = 0.625 and resolution (30,15,15), an 
apparently chaotic edge trajectory was followed for 0{hOQD jU) during which 5 different 
phases of relative recurrence - labelled A2, i?2, C2, D2 and E2 in chronological order - 
were detected. Interestingly, this edge trajectory, although seemingly chaotic, is much 
smoother than the corresponding (same Re and a) edge trajectory for mp — 1: compare 
Figures [T] and [TOl Of the 5 initial conditions for a Newton-Krylov search, r^in only 
converged to O(10~") for D2 and E2 {r^in for A2,B2 and C2 stagnated at 0(10-^)). 

The exact recurrence states found, i?2_0.625 and £'2_0.625 (see Figure[TT|), correspond 
to two new types of travelling wave branches distinct from those originally reported in 
Faisst & Eckhardt (2003) and Wedin & Kerswell (2004). i:»2_0.625 is especially interesting 
as it is the first TW found in pipe flow which is not S-symmetric. It does, however, possess 
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a mirror symmetry 

: (m, V, w,p){s, (j), z) {u, -V, w,p){s, 2ip - 0, z) (3.5) 

where ip « 45° (or equivalently « —45° because of the imposed R2-symmetry) for the 
snapshot shown in Figure [TT] In contrast, i?2_0.625 has three symmetries: two unimposed 
- S and Z±7r/4 (where (f) — defines the plane of shift-&-reflect symmetry) - and one 
imposed - R2. i?2_0.625 is actuahy one member of a whole new class of /lig/iZj/ -symmetric 
TWs which have subsequently been uncovered (the TV-class, see Pringle et al. 2008). These 
new waves are significant because they appear much earlier in Re than the originally 
discovered TWs (Faisst & Eckhardt 2003, Wedin & Kerswell 2004) and their upper 
branch solutions look more highly nonlinear. 

3.6. Local relative attractors within the 112 -symmetric subspace 

In this subsection, the flow is again restricted to be R2-symmetric but the pipe is halved 
in length (a = 1.25 so L « 2.5 D) and Re reduced to 2400 (resolution is (50,20,8)). At 
these settings, Kerswell & Tutty (2007) identified a lower branch TW - 25_1.25 in their 
nomenclature - which was on the edge with only one unstable direction normal to it. 
This TW is therefore a local attractor for trajectories confined to the edge. A leading 
question is whether it is also a global attractor on the edge. 

A starting point was randomly chosen along a turbulent trajectory, and an edge tra- 
jectory generated for 500 D/U using the method described above. The energy signal is 
displayed in Figure [T^] and the corresponding recurrence signal rmin (t) is displayed in 
Figure fT3l Two successive dips in rmin can be seen at times 50 D/U and 90 D/U later 
corresponding to inexact recurrences labelled respectively A3 and B3. When used to 
initialise a Newton-Krylov search, both led to a value of 0(10^^^) for rmm. The two 
converged exact states correspond to one unique travelling wave solution, differing only 
by a small shift in the azimuthal direction of less than 4 degrees. As a result, we label 
only the first A3A.25: see Figure [TTl The TW A3A.25 has all the same symmetries as 
i?2_0.625 but possesses a distinctly different structure consisting of 4 low-speed streaks 
near the walls separated by 4 high-speed streaks. Four other high-speed streaks of lesser 
amplitude and smaller size, and with a much more steady structure, can be found closer 
to the pipe axis. This new TW thus has a richer structure than the R2 branch of so- 
lutions found in Faisst & Eckhardt (2003) and Wedin & KersweU (2004). Just as for 
£'2_0.625, subsequent investigations have revealed that A3A.25 is but one member of 
another highly-symmetric class of TWs (the M-class, Pringle et al. 2008) characterised 
by this distinctive double-layer structure. 

After a time 200 D/U, the energy of the edge trajectory reaches a constant plateau 
and the value of rmin decreases exponentially to O(10~'*) by 300 D/U. This indicates 
that the trajectory is converging to an end state. Newton-Krylov was used to accelerate 
the convergence of the trajectory by taking the state at 300 D/U and reducing rmin to 
0(10""'^^). The spatial structure of the initial and final state of this procedure, however, is 
indistinguishable. The final state is another travelling wave solution, C3_1.25 (see Figure 
[TT|) . distinct from 2&_1.25 of KersweU & Tutty (2007). Stability analysis of C3.1.25 reveals 
only one unstable direction which has to be normal to the edge (see Pringle et al 2008 
for details). Hence C3_1.25 is an attractor in the edge, and both C3_1.25 and 26_1.25 can 
only be local attractors there. Closer examination of i?2_0.625 and C3_1.25, which look 
very similar, reveals that they are members of the same TW family distinguished only 
by their different wavelengths. 
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3.7. Travelling waves on the Edge 

There is no guarantee that starting guesses taken from the edge for the Newton-Kryfov 
method wiU converge to states also on the edge. The calculations of subsection 3.4 have, 
however, demonstrated this for Al_0.625 (the asymmetric TW) thereby confirming earlier 
speculation (Schneider et al 2007b, Pringle & Kerswell 2007). To explicitly check this for 
the other TWs found here, each was taken as an initial condition with various very small 
perturbations added to start a series of runs. In all cases, two contrasting behaviours were 
evident: one perturbation could be found which led to relaminarisation while another led 
to a turbulent episode: see Figure [T31 This confirms that all the TWs found here are on 
the laminar/turbulent boundary for the particular parameter settings tested (shown in 
Table 1). 



4. Discussion 

A summary of the main findings of this investigation are listed below. 

• The laminar-turbulent boundary S in a short pipe seems to have a chaotic attractor 
as observed by Schneider et al. (2007b). 

• Although S has many unstable lower branch TWs embedded within it (e.g. Kerswell 
& Tutty 2007), only one - the asymmetric TW of Pringle & Kerswell (2007) (and rotated 
derivatives of it) - was approached by edge trajectories studied here. 

• A new asymmetric TW which has a small but finite azimuthal phase speed has been 
discovered. 

• Some evidence was gathered for a "relative" homoclinic connection between the 
asymmetric TW and the same wave with a different orientation. 

• Calculations within a R2-symmetric subspace have revealed new branches of R2- 
symmetric TW solutions. One (Z32_0.625) is the first TW found not to be shift-&-reflect 
symmetric in pipe flow, while the others are members of new classes of /lig/iZj/ -symmetric 
TWs (Pringle et al. 2008). 

• The laminar-turbulent boundary restricted to R2-symmetry has at least two simple 
local attractors in the form of TWs at Re = 2400. 

Our numerical experiments have shed some light on the dynamical structure of E 
in pipe flow. When no symmetry is imposed, and the pipe is approximately 5 D long 
(a = 0.625), trajectories which neither relaminarise nor become turbulent all look chaotic 
and visit some exact recurrent states in a transient manner. In all cases studied, these ex- 
act solutions appear to be speciflcally the asymmetric TWs found by Pringle & Kerswell 
(2007) even though many other lower branch TWs are embedded in S (Kerswell & Tutty 
2007). No RPOs were found even though the numerical method was sufficiently general 
to find them, in addition to TWs. Whether this is because the RPOs are rarely visited 
or just more difficult to isolate is unclear. It is tempting to conclude that chaos on the 
edge of pipe flow is structured around a set of unstable saddle points (the TWs) linked 
together by heteroclinic or sometimes relative heteroclinic connections. An approxima- 
tion to such a connection has been shown in subsection 13.41 The resulting dynamical 
structure is likely to be a heteroclinic (or homoclinic) tangle, an efficient mechanism to 
produce chaotic trajectories in phase-space with the set of TWs acting as a skeleton of 
the tangle. Seen from this point of view, transition to turbulence from a given initial con- 
dition depends on the position of the initial state, in phase space, relative to the stable 
manifold of each exact state embedded in E. When a trajectory approaches a travelling 
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TW 


mo 




C 




status 


yll_0.625 


1 


2875 


1.55494 





known 


A'_0.625 


1 


2875 


1.51956 


-1.704.10-^ 


new 


D2_0.625 


2 


2875 


1.53382 





new 


^3_1.25 


2 


2400 


1.23818 





new 


C3_1.25 


2 


2400 


1.55064 





new 



Table 1. Parameters, axial propagation speed C (in units of U) and azimuthal propagation 
speed Ce in radU/D and status of all travelling wave solutions found in this paper (i52_0.625 is 
in the same TW family as C3_1.25 and hence not listed separately). 



wave solution belonging to S, its relative position to the boundary determines on which 
side of the edge the trajectory escapes, resulting in either relaminarisation or a turbulent 
transient (see sketch in Figure [5]). 

When the flow is artificially restricted to be R2-symnietric at Re = 2875, the trajec- 
tory remains chaotic wandering sufficiently close to distinct TW states (£'2_0.625 and 
£'2_0.625) to be captured by the Newton-Krylov scheme. By halving the pipe length to 
L « 2.5 D and reducing Re to 2400, multiple local attractors appear in E. A randomly- 
started trajectory, after a chaotic transient and approaches to states ^3_1.25 and 53-1.25, 
is attracted towards the TW solution C3_1.25. The fundamental difference between 
C3_1.25 and all the other states mentioned so far lies in the number of unstable eigenval- 
ues. When a TW has two or more unstable eigenvalues in a given subspace, the dimension 
of the intersection between its unstable manifold and the hypersurface E is reduced by 
one but remains at least one. Hence such a state remains a saddle point on E: edge trajec- 
tories enter its vicinity along its stable manifold and escape along the unstable manifold. 
In contrast, when a state has only one unstable direction, this is necessarily normal to 
E, so that the state becomes an attractor for dynamics restricted to the edge (a relative 
attractor). We have shown that C3_1.25 is one such TW, but there is at least one other, 
26_1.25 of Kerswell & Tutty (2007). Therefore we expect both these two states to be 
only local attractors on E rather than global ones. In the cases of plane Poiseuille and 
Couette flow, a lower-branch solution exists which has only one unstable direction and 
is therefore a relative attractor (Toh & Itano 1999, Itano & Toh 2001, Wang et. al. 2007, 
Schneider et al., 2008, Viswanath 2008). Pipe flow is different as TWs all seem to have at 
least two unstable directions when there is no imposed discrete rotational symmetry. For 
example, the 26_1.25 TW has 1 harmonic (R2-symmetric) unstable eigenfunction and 
one subharmonic (only Ri-symmetric) unstable eigenfunction. 

While discrete rotational symmetry constraints might appear artificial, they have nev- 
ertheless served as a useful device for discovering new exact recurrent solutions (see Table 
1). The method developed here based on edge tracking, recurrence analysis and use of a 
Newton-Krylov algorithm, is a general approach for finding new exact recurrent solutions 
in any flow situation possessing subcritical behaviour. Importantly, the present method 
naturally selects the states that are most likely to be visited and does not presuppose 
anything about their spatial structure (modulo the symmetries imposed on the flow). The 
use of the scalar function rmim coupled with a Newton-Krylov solver, can also be used 
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to search for more complex solutions such as relative periodic orbits, whether located 
on the laminar-turbulent boundary or embedded in fully developed turbulence. This is 
currently underway. 

May open questions remain regarding E in pipe flow. Is S really a hypersurface or can 
it have a more fractal structure? Why does the flow keep approaching the asymmetric 
TW or rotations of it, and not any of the other lower branch TWs known to exist within 
the same parameters range? Finally, what are the properties of S for longer pipes, where 
localised turbulent 'puffs' structures exist? Preliminary work for the extended pipe in a 
reduced model has already revealed an interesting localised structure as the attracting 
state (WilUs & Kerswell 2008b). 



Acknowledgements: 

Wc would like to thank C. Pringlc for helping with the velocity profiles of the TW so- 
lutions and B. Eckhardt for stimulating discussions about this topic. Y.D. was supported 
by a Marie-Curie Intra-European Fellowship (grant number MEIF-CT-2006-024627) and 
A.W. by the EPSRC (grant number GR/S76144/01). 



REFERENCES 

Brown, P. N. & Saad, Y. 1990 Hybrid Krylov methods for nonlinear systems of equations 

Siam J. Sci. Stat. C'omput, 11, 450-481 
Dennis, J. E. Jr. & Schnabel, R. E. 1995 Numerical Methods for Unconstrained Optimization 

and Nonlinear Equations, SIAM 
EiSENSTAT, C. & Walker, H. 1996 Choosing the forcing terms in an inexact Newton method 

SIAM J. Sci. Comput, 17, 16-32 
Faisst, H. & ECKIIARDT, B. 2003 Travelling waves in pipe flow Phys. Rev. Letters, 91,224502 
Faisst, H. & ECKIIARDT, B. 2004 Sensitive dependence on initial conditions in transition to 

turbulence J. Fluid Mech., 504, 343-352 
Gibson, J. F., Halcrow, J., Cvitanovic, P. 2008 The geometry of state-space in plane 

Couette flow, submitted to J. Fluid Mech. 
Hagen, G. H. L. 1839 Uber die Bewegung des Wassers in engen zylindrischen Rohren Poggen- 

dorfs Annalen der Physik und Chemie, 16, 423 
Hamilton, J.M., Kim, J. & Waleffe, F. 1995 Regeneration mechanisms of near-wall turbu- 
lence structures. J. Fluid Mech, 287, 317-348 
Hop, B., van Doorne, C. W. H., Westerweel, J., Nieuwstadt, F. T. M., Faisst, H., 

Eckhardt, B., Wedin, H., Kerswell, R.R. & Waleffe, F. 2004 Science, 305, 1594- 

1597 

HoF, B., VAN Doorne, C. W. H., Westerweel, J. & Nieuwstadt, F. T. M. 2005 Phys. 

Rev. Lett, 95, 214502 

Itano, T. & TOH, S. 2001 The dynamics of bursting process in wall turbulence J. Phys. Soc. 
Japan 70, 703-716 

Jimenez, J. & Moin P. 1991 The minimal flow unit in near- wall turbulence, J. Fluid Mech., 
225, 213-240 

Joseph, D. D. & Carmi, S. 1969 Stability of Poiseuille flow in pipes, annuli, and channels 

Quart. Appl. Math., 26, 575-599 
Kawahara, G. 2005 Laminarization of minimal plane Couottc flow: Going beyond the basin of 

attraction of turbulence Phys. Fluids, 17, 041702 
Kerswell, R. R. 2005 Recent progress in understanding the transition to turbulence in a pipe 

Nonlinearity, 18, R17-R44 



Transition in pipe flow 



13 



Kerswell, R. R. & TuTTY, O.R. 2007 Recurrence of Travelling Waves in Transitional Pipe 
Flow J. Fluid Mech., 584, 69-102 

Marques, F. 1990 On boundary conditions for velocity potentials in confined flows; Application 
to Couette flow Phys. Fluids, A2 , 729-737 

Mellibovsky, F. & Meseguer, a. 2007 Pipe flow dynamics on the critical threshold Proceed- 
ings of the 15th. Int. Couette- Taylor Workshop (ed. Mutabazi, I.), Le Havre, France 

Nagata, M. 1990 Three dimensional finite amplitude solutions in plane Couette flow: bifurcation 
from infinity. J. Fluid Mech., 217 519-527 

NiKlTiN N. 2007 Spatial periodicity of spatially evolving flows caused by inflow boundary con- 
dition, Phys. Fluids 19, 091703 

Peixinho, J. & MULLIN, T. 2006 Decay of turbulence in pipe flow Phys. Rev. Lett., 96, 094501 

Peixinho, J. & MuLLiN, T. 2007 Finite-amplitude thresholds for transition in pipe flow J. 
Fluid Mech., 582, 169-178 

Pfenniger, W. 1961 Transition in the inlet length of tubes at high Reynolds numbers Boundary 
Layer and Flow Control (ed. GV Lachman), 970 

Poiseuille, J. L. M. 1840 Recherches experimentales sur le mouvement des liquides dans les 
tubes de tres petits diametres CR. Acad. Sci., 11, 961 

Pringle, C. C. T. & Kerswell, R. R. 2007 Asymmetric, helical and mirror-symmetric trav- 
elling waves in pipe flow Phys. Rev. Letters, 99, 074502 

Pringle, C. C. T., Duguet, Y. & Kerswell, R. R. 2008 Highly-symmetric travelling waves 
in pipe flow Phil Trans Roy Soc, submitted (http://arxiv.org/abs/0804.4854) 

Priymak, V. G., MiYAZAKi, T 2004 Direct numerical simulation of equihbrium spatialUy- 
localized structures in pipe flow Phys. Fluids, 16, 4221-4234 

Reynolds, O. 1883 An experimental investigation of the circumstances which determine 
whether the motion of water shall be direct or sinuous and of the law of resistance in 
parallel channels Phil. Trans. Roy. Soc, 174, 935-982 

Saad, Y. & SCHULTZ, M.H. 1986 GMRES: a generalized minimal residual method for solving 
nonsymmetric linear systems Stam J. Sci. Stat. Comput, 7, 856-869 

Schneider, T. M., Eckhardt, B. & Vollmer, J. A. 2007 Statistical analysis of coherent 
structures in transitional pipe flow PRE, 75, 066313 

Schneider, T. M., Eckhardt, B. & Yorke, J. A. 2007 Turbulence transition and the edge 
of chaos in pipe flow Phys. Rev. Lett., 99, 034502 

Schneider, T. M., Gibson J. F., Lagha, M., De Lillo, F. & Eckhardt, B. 2008 Laminar- 
turbulent boundary in plane Couette flow preprint (http://arxiv.org/abs/0805.1015) 

Skufca, J.D., Yorke, J. A. & Eckhardt, B. 2006 Edge of Chaos in a Parallel Shear Flow 
Phys. Rev. Letters, 96, 174101 

Toh, S. & Itano, T. 1999 Low-dimensional dynamics embedded in a plane Poiseuille flow 
turbulence: Traveling-wave solution is a saddle point? Proc. lUTAM Symp. on Geometry 
and Statistics of Turbulence (ed. Kambe, T.), Kluwer 

Viswanath, D. 2007 Recurrent motions within plane Couette turbulence J. Fluid Mech., 580, 
339-358 

Viswanath, D. 2008 The dynamics of transition to turbulence in plane Couette Mathematics 

and Computation - A Contemporary view. Springer, 3, Abel Symposia 
Waleffe, F. 1997 On the self-sustaining process in shear flows Phys. Fluids, 9 883-900 
Waleffe, F. 1998 Three dimensional coherent states in plane shear flows Phys. Rev. Lett., 81 
4140-4143 

Waleffe, F. 2001 Exact Coherent Structures in channel flow J. Fluid Mech., 435 93-102 
Waleffe, F. 2003 Homotopy of exact coherent structures in plane shear flows Phys. Fluids, 15 
1517-1534 

Wang, J., Gibson, J., Waleffe, F. 2007 Lower branch coherent states in shear flows : tran- 
sition and control Phys. Rev. Lett. 98, 204501. 

Wedin, H. & Kerswell, R. R. 2004 Exact coherent structures in pipe flow: travelling wave 
solutions J. Fluid Mech., 508, 333-371 

Willis, A. P. & Kerswell, R. R. 2007 Critical behaviour in the relaminarisation of localised 
turbulence in pipe flow Phys. Rev. Lett. 98, 14501. 

Willis, A. P. & Kerswell, R. R. 2008a Coherent structures in local and global pipe turbu- 
lence Phys. Rev. Lett. 100, 124501 



14 



Y. Duguet, A. P. Willis and R.R. Ker swell 



0.01 



0.001 



le-04 



1e-05 




200 250 
time 



450 



Figure 1. Energy contained in the axially-dependent modes for [a, Re, mo) = (0.625,2875,1). 
The thick hne indicates the edge trajectory and the thinner lines nearby trajectories which either 
relaminarise (energy decreases) or becomes turbulent (energy increases to a higher level). Time 
is in units of D/U . 
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Figure 2. Axial perturbation velocity signal at two different locations in the pipe: 
(s = 0.66, 6* = 0,2 = 0) and (s = 0.76, 6* = 0,z = 0.5 D). The signal is taken from the edge 
trajectory in Section [3.11 Dotted arrows indicate where the velocity signal shows approximate 
recurrences. 

1.2 I ^ ^ ^ 1 1 




50 100 150 200 250 

time 



Figure 3. Typical profiles of the residual function r{t) starting from snapshots of the edge 
trajectory of Section WT\ The subscript i has been suppressed (time in units of D/U). 
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LAMINAR FLOW 



Figure 5. Schematic view of phase-space. The surface E separates initial conditions which 
relaminarise from those which become turbulent. An edge trajectory visiting three states A, B 
and C is shown schematically. The dynamics on the manifolds transverse to S are shown by 
trajectory diverging towards either the laminar state or the turbulent state. 
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Figure 6. Starting guesses for the Newton-Krylov algorithm (top left Al, top right Bl) and 
converged states (bottom left yll_0.625, bottom right _B1_0.625) on the edge trajectory described 
in i)3.2l Each subfigure represents a snapshot across the pipe. Contours indicate the axial velocity 
difference from the underlying laminar flow (light/dark indicating faster /slower moving fluid) 
and the arrows represent the cross-stream velocity (length proportional to speed). Maximum 
norm of the (u, v) cross- velocity is 0.0128C/ and the ax;ial velocity differential w is in the range 
±0.17?7. The TWs j41_0.625 and i?l_0.625 are exactly the same state modulo an azimuthal shift. 
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Figure 7. A sketch of the heterocUnic connnection joining two saddle points A and B. ei and 
62 are the two unstable eigendirections of A, used for shooting by varying the shooting angle, 
in order to optimise the approach to B. 
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Figure 8. Energy contained in the axially dependent modes for Re = 2875, mo = 1, a = 0.625. 
The thick hne indicates the edge trajectory and the thinner lines nearby trajectories which either 
relaminarise (energy decreases) or becomes turbulent (energy increases to a higher level). Time 
is in units oi D/U 
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Figure 9. Recurrence signal for the heteroclinic connection of H3. 41 versus time (in D/U). 

{a,Re,mo) = (0.625,2875,1). 
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Figure 10. Energy contained in the axially-dependent modes on the edge for 
(a, Re, mo) = (0.625, 2875, 2). The thick hne indicates the edge trajectory and the thinner hnes 
nearby trajectories which cither rclarninarisc (energy decreases) or becomes turbulent (energy 
increases to a higher level). Time is in units oi D/U 
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Figure 11. Converged states D2_0.625 (top, left), £;2_0.625 (top, right), yl3_1.25 (bottom, left) 
and C3_1.25 (bottom,riglit). Contours indicate the axial velocity difference from the underlying 
laminar flow (light/dark indicating faster/slower moving fluid) and the arrows represent the 
cross-stream velocity (length proportional to speed). Across the four snapshots, the maximum 
cross-stream speed is 0.0143(7, and the axial velocity perturbation w is in the range ±0.17(7. 
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Figure 12. Energy contained in the axially-dependent modes for {a, Re, mo) = (1.25,2400,2). 
The thick line indicates the edge trajectory and the thinner lines nearby trajectories which either 
relaminarise (energy decreases) or becomes turbulent (energy increases to a higher level). Time 
is in units oi D/U 



0.1 



0.01 



A3 



B3 



0.001 



1e-04 




350 



time 



Figure 13. Recurrence signal for the attracting state of Section [3. 51 versus time expressed in 
D/U units, {a, Re, mo) = (1.25,2400,2). 
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Figure 14. Diverging trajectories demonstrating tlie 'edge credentials' of the new TWs found 
here. Two typical trajectories starting close to each TW axe shown with one decaying down 
smoothly to the laminar state (Esa = 0) and another increasing up to the turbulent state. D2 
represents D2_0.625, A3 is A3A.25, A is yll_0.625 and C3 is C3_1.25 with ± used to distinguish 
the starting perturbation. 



